Table of Contents
- 1 2-1
- 2 2-2
- 3 2-3 (Monthly Unemployment Rate)
- 4 2-4 (CRSP Decile Index)
- 5 2-5 (IBM)
- 6 2-6 (Monthly Power Consumption)
- 7 2-7 (Equal-Weighted Index)
- 8 2-8 (S&P)
- 9 2-9 (IBM)
- 10 2-10 (Aaa and Baa Bond Yields)
- 11 2-11 (Aaa and Baa Bond Yields)
- 12 2-12 (Aaa and Baa Bond Yields)
- 13 2-13 (CRSP Equal-Weighted Index)
- 14 2-14 (S&P 500 Spot and Futures)
- 15 2-15 (Quarterly GDP Implicit Prce Deflator)
- 16 End
In [138]:
print(R.version)
file.path(R.home("bin"), "R")
In [1]:
library(IRdisplay) # display
library(fUnitRoots) # adfTest
library(TSA) # eacf
In [511]:
# rm(plot_arima_forecast_fig)
In [631]:
# detach("package:AFTSCode", unload = TRUE) # First detach the package
# unloadNamespace("AFTSCode")
In [632]:
library(AFTSCode)
2-1¶
- From Sec. 2.5.4
- $E[R_{101}|F_{100}] = c_0 - \theta_1 a_{100} = 0.2*0.01 = 0.002$
- $Var[R_{101}|F_{100}] = \sigma_{a}^2 = 0.025^2 = 0.000625$
- $E[R_{102}|F_{100}] = c_0 = 0$
- $Var[R_{102}|F_{100}] = (1+\theta_1^2)\sigma_{a}^2 = (1+0.2^2)*0.025^2 = 0.00065$
- From Sec. 2.5.1
- $\rho_1 = -(-0.2)/(1+0.2^2) = 0.192307692307692$
- $\rho_2 = 0$
2-2¶
- From Sec. 2.4.1
- $E[r_t]=\frac{\phi_0}{1-\phi_1}=0.01/(1-0.2)=0.0125$
- $Var[r_t]=\frac{\sigma_a^2}{1-\phi_1^2}=0.02/(1-0.2^2)=0.0208333333333333$
- From Sec. 2.4.2
- $\rho_0=1$
- $\rho_1=\phi_1\rho_0=0.2$
- $\rho_2=\phi_1\rho_1=0.04$
- From Sec. 2.4.4
- $E[r_{101}|F_{100}]=\phi_0+\phi_1*r_{100}=0.01+0.2*(-0.01)=0.008$
- $Var[r_{101}|F_{100}]=\sigma_a^2=0.02$
- $E[r_{102}|F_{100}]=\phi_0+\phi_1E[r_{101}|F_{100}]=0.01+0.2*0.008=0.0116$
- $Var[r_{102}|F_{100}]=Var[\phi_1 r_{101}+a_{102}|F_{100}]=Var[\phi_1 (\phi_0+\phi_1 r_{100}+a_{101})+a_{102}|F_{100}]=(1+\phi_1^2)\sigma_a^2=(1+0.2^2)*0.02=0.0208$
2-3 (Monthly Unemployment Rate)¶
- Interesting observations.
- The
diff(log(unem_rate))shows first strong momentum, second strong mean-reverting, and serial correlation inpacf. Seeplot_pacf_acfresults.
- The
In [554]:
da = read.table("../AFTS_sol/data/m-unrate.txt", header = TRUE)
da[1:5,]
In [555]:
unem_rate = da$Rate
lg_unem_rate = log(unem_rate)
unem_rate_ts = ts(unem_rate, frequency = 12, start = c(1948,1))
lg_unem_rate_ts = ts(lg_unem_rate, frequency = 12, start = c(1948,1))
In [143]:
plot_time_fig(unem_rate_ts, main = "Unemployment Rate", xlab = "Year")
pacf(unem_rate)
In [22]:
plot_time_fig(lg_unem_rate_ts, main = "log(Unemployment Rate)", xlab = "Year")
pacf(lg_unem_rate)
In [557]:
plot_pacf_acf(lg_unem_rate, freq = 12, lag.max = 60)
Unit-root test¶
We can reject the null-hypothesis from ADF test. In other words, the process doesn't have unit-root.- Seems like we cannot reject he null hypothesis when we use large enough lags in the ADF testing.
- Also from the ACF plots of the unemployment rate before and after taking difference, we can see that we need to take difference to cancel the strong serial correlation.
- From the ACF plot of the unemployment rate after taking difference, we see some damping sine and cosine wave, indicating that there could be some business cycle (pp42).
In [34]:
# "lags=12" is chosen from PACF plot
adf_test_res = adfTest(unem_rate_ts, lags = 12, type = c("ct"))
adf_test_res
In [31]:
lag.max = 60
plot_acf(da = unem_rate, lag.max = lag.max, main = "Unemployment Rate")
pacf(unem_rate, lag.max = lag.max)
In [ ]:
In [35]:
diff_unem_rate = diff(unem_rate)
diff_unem_rate_ts = ts(diff_unem_rate, frequency = 12, start = c(1948,2))
In [36]:
plot_acf(da = diff_unem_rate, lag.max = lag.max, main = "diff(Unemployment Rate)")
pacf(diff_unem_rate, lag.max = lag.max)
In [144]:
plot_acf(da = diff(unem_rate, 12), lag.max = lag.max, main = "diff(Unemployment Rate, 12)")
pacf(diff(unem_rate, 12), lag.max = lag.max)
Order determination¶
- Try to use ARIMA(4, 1, 5) model.
In [172]:
perform_and_print_eacf <- function(da, ar.max, ma.max) {
eacf_obj <- eacf(da, ar.max = ar.max, ma.max = ma.max)
eacf_stats_tb = format(as.data.frame(eacf_obj$eacf), digits = 3)
names(eacf_stats_tb) <- seq(from = 0, to = ma.max)
display(eacf_stats_tb)
display(eacf_obj$symbol)
# pp67, asymptotic standard error of EACF
display(2/sqrt(length(da)))
c(eacf_obj, eacf_stats_tb)
}
In [173]:
unem_rate_eacf_res <- perform_and_print_eacf(diff_unem_rate, ar.max = 25, ma.max = 12)
Fit an ARIMA(3,1,5) model¶
In [108]:
unem_rate_mod = stats::arima(unem_rate_ts, order = c(3,1,5), include.mean = F)
unem_rate_mod
In [110]:
unem_rate_mod_1 = arima(unem_rate_ts, order = c(3,1,5))
unem_rate_mod_1
In [111]:
diff_unem_rate_mod = arima(diff_unem_rate_ts, order = c(3,0,5), include.mean = F)
diff_unem_rate_mod
In [112]:
plot_time_fig(unem_rate_mod$residuals, main = "Residuals")
plot_acf(unem_rate_mod$residuals, main = "ACF of Residuals")
Box-Ljung test¶
- pp33 very briefly talks about how to choose lag in the
Box.test - From the following testing result, we cannot reject the null-hypothesis. So the model specified is adequate.
In [113]:
length(unem_rate); log(length(unem_rate))
In [114]:
Box.test(unem_rate_mod$residuals, lag = 7, type = 'Ljung')
In [115]:
Box.test(unem_rate_mod$residuals, lag = 12, type = 'Ljung')
- Diagnose time-series ARIMA model
In [161]:
help(tsdiag)
In [162]:
par(bg = 'white')
tsdiag(unem_rate_mod, gof.lag=36)
Check business cycle¶
- pp42
- Compute the average length of business cycles
- Doesn't matter if using the characteristic roots or the inverse of them.
- Unit is "month". Average business cycle is about 14.2 months, which seems to be shoter than expected from the time plot above.
- The business cycle is very sensitive to the order chosen in fitting ARIMA model. When I use ARIMA(2,0,5) or ARIMA(4,0,5), the business cycle doens't exist or too short.
In [69]:
find('arima')
In [79]:
help("arima", package = "TSA")
In [80]:
help("arima", package = "stats")
In [116]:
sqrt(unem_rate_mod$sigma2)
In [117]:
unem_rate_mod$coef
In [118]:
ar_poly = c(1, -unem_rate_mod$coef[1:3]) # Characteristic equation for AR process
roots = polyroot(ar_poly)
roots
In [119]:
for (i in 1:3) {
print(paste("====", i))
print(as.numeric(Mod(1-unem_rate_mod$coef[1]*roots[i]-unem_rate_mod$coef[2]*roots[i]^2-unem_rate_mod$coef[3]*roots[i]^3-unem_rate_mod$coef[4]*roots[i]^4)))
print(as.numeric(Mod(roots[i])))
}
- Modulus of roots
In [120]:
roots; Mod(roots)
In [121]:
1/roots; Mod(1/roots)
- Compute the average length of business cycles
- Doesn't matter if using the characteristic roots or the inverse of them.
In [122]:
roots[1]; Re(roots[1]); Im(roots[1])
In [125]:
k = 2*pi/acos(Re(roots[1])/Mod(roots[1]))
k
In [126]:
k = 2*pi/acos(Re(1/roots[1])/Mod(1/roots[1]))
k
Forecast¶
In [128]:
length(unem_rate_ts); unem_rate_ts
In [140]:
npts = 8
eotr = length(unem_rate_ts)
h = 4
freq = 12
order = c(3,1,5)
fixed = NULL
seasonal = list(order = c(0L, 0L, 0L), period = NA)
unem_rate_fc_res = plot_arima_forecast_fig(
da_ts=unem_rate_ts, eotr=eotr, h=h, npts=npts, frequency=freq,
order=order, seasonal=seasonal, fixed=fixed, method='ML',
include.mean=F, transform.pars=TRUE,
main="Forecasts from ARIMA(3,1,5) for unemp_rate",
xlab="Year", ylab="Unemp-Rate", ylim=c(5, 11)
)
Alternative solutions (only use AR model)¶
In [147]:
par(bg = 'white')
# decay slowly
acf(unem_rate)
par(mfrow = c(2, 1))
acf(diff(unem_rate))
pacf(diff(unem_rate))
acf(diff(unem_rate, 12))
pacf(diff(unem_rate, 12))
cat("order =", ar(unem_rate, method = 'mle')$order, "\n")
m1 = arima(unem_rate_ts, order = c(11, 0, 0))
m1
# Follow the pp51 to modify and refine the model based on previous fit
m1 = arima(unem_rate_ts, order = c(11, 0, 0), fixed = c(NA, NA, 0, 0, 0, NA, 0, 0, 0, NA, NA, NA))
m1
# FIXME
m2 = arima(unem_rate_ts, order = c(2, 1, 1), seasonal = list(order = c(1, 0, 1), period = 12))
m2
tsdiag(m1, gof = 36)
tsdiag(m2, gof = 36)
predict(m1, 4)
predict(m2, 4)
In [149]:
ar_poly_1 = c(1, -m1$coef[1:11]) # Characteristic equation for AR process
roots_1 = polyroot(ar_poly_1)
roots_1
In [155]:
Im(roots_1[8])
In [156]:
for (i in 1:length(roots_1)) {
if (abs(Im(roots_1[i]))>1e-8) {
print(i)
print(roots_1[i])
print(2*pi/acos(Re(roots_1[i])/Mod(roots_1[i])))
}
}
In [151]:
ar_poly_2 = c(1, -m2$coef[1:2]) # Characteristic equation for AR process
roots_2 = polyroot(ar_poly_2)
roots_2
2-4 (CRSP Decile Index)¶
- Data are monthly simple returns.
- Interesting observations.
- The
log(dec2_rtn)shows medium momentum correlation inpacfandacfat lag-1 and lag-12. Thelog(dec10_rtn)shows almost no serial correlation. Seeplot_pacf_acfresults.
- The
In [546]:
da = read.table("../AFTS_data/Ch02/m-deciles08.txt", header = TRUE)
dim(da); da[1:5,]
In [547]:
dec2 = da$CAP2RET
dec2_ts = ts(da$CAP2RET, frequency = 12, start = c(1970, 1))
dec10 = da$CAP10RET
dec10_ts = ts(da$CAP10RET, frequency = 12, start = c(1970, 1))
In [5]:
plot_time_fig(dec2_ts, main = "CRSP Decile Index", xlab = "Year")
lines(dec10_ts, type = , col = 'red', lty = 2)
legend(
"topright",
legend = c("Decile 2", "Decile 10"),
col = c("black", "red"),
lty = c(1, 2),
pch = c(8, NA)
)
In [548]:
plot_pacf_acf(dec2, freq = 12, lag.max = 36)
In [551]:
plot_pacf_acf(log(1+dec2), freq = 12, lag.max = 36)
In [549]:
plot_pacf_acf(dec10, freq = 12, lag.max = 36)
In [552]:
plot_pacf_acf(log(1+dec10), freq = 12, lag.max = 36)
(a) Test autocorrelations¶
- Reject the null-hypothesis (i.e., the first 12 lags have no serial correlation) for Decile 2. So for Decile 2, the first 12 lags have serial correlation.
- Cannot reject the null-hypothesis for Decile 10. So for Decile 10, the first 12 lags have no significant serial correlation.
In [6]:
par(mfrow = c(2, 1), bg = 'white')
acf(da$CAP2RET)
pacf(da$CAP2RET)
Box.test(da$CAP2RET, lag = 12, type = 'Ljung')
In [7]:
par(mfrow = c(2, 1), bg = 'white')
acf(da$CAP10RET)
pacf(da$CAP10RET)
Box.test(da$CAP10RET, lag = 12, type = 'Ljung')
(b) Fit a model and perform model checking¶
In [8]:
dec2_eacf_res <- perform_and_print_eacf(da$CAP2RET, ar.max = 25, ma.max = 12)
- Fit
ARIMA(1,0,1)model, b/c only this model has all coefs significant.
In [10]:
dec2_mod_1 <- arima(da$CAP2RET, order = c(1, 0, 1), include.mean = T)
dec2_mod_1; Box.test(dec2_mod_1$residuals, lag = 12, type = 'Ljung');
par(bg = 'white')
tsdiag(dec2_mod_1, gof.lag = 36)
In [11]:
# Follow the pp51 to modify and refine the model based on previous fit
dec2_mod_2 <- arima(da$CAP2RET, order = c(0, 0, 1), include.mean = T)
dec2_mod_2; Box.test(dec2_mod_2$residuals, lag = 10, type = 'Ljung');
par(bg = 'white')
tsdiag(dec2_mod_2, gof.lag = 36)
In [12]:
dec2_mod_3 <- arima(da$CAP2RET, order = c(4, 0, 3), include.mean = T)
dec2_mod_3; Box.test(dec2_mod_3$residuals, lag = 12, type = 'Ljung');
par(bg = 'white')
tsdiag(dec2_mod_3, gof.lag = 36)
In [13]:
dec2_mod_3 <- arima(
da$CAP2RET,
order = c(4, 0, 3),
# Follow the pp51 to modify and refine the model based on previous fit
fixed = c(0, NA, 0, NA, 0, NA, 0, NA),
include.mean = T
)
dec2_mod_3; Box.test(dec2_mod_3$residuals, lag = 12, type = 'Ljung');
par(bg = 'white')
tsdiag(dec2_mod_3, gof.lag = 36)
(c) Forecast using ARMA model¶
In [22]:
npts = 12
eotr = length(dec2_ts)-npts
h = npts
freq = 12
order = c(0,0,1)
fixed = NULL
seasonal = NULL
dec2_fc_res = plot_arima_forecast_fig(
da_ts=dec2_ts, eotr=eotr, h=h, npts=npts, frequency=freq,
order=order, seasonal=seasonal, fixed=fixed, method='ML',
include.mean=T, transform.pars=TRUE,
main="Forecasts from ARIMA(0,0,1) for Decile 2 index",
xlab="Year", ylab="Unemp-Rate"# , ylim=c(5, 11)
)
dec2_fc_tb = comb_forecast_res(dec2_fc_res, dec2_ts, eotr, freq)
dec2_fc_tb
In [17]:
class(dec2_fc_res); attributes(dec2_fc_res); dec2_fc_res$model
In [20]:
help("predict")
In [21]:
predict(dec2_fc_res$model, 12)
In [23]:
dec2_fc_tb
2-5 (IBM)¶
- Data are daily simple returns.
- There is some long-range dependency in daily IBM simple returns.
In [561]:
da = read.table("../AFTS_sol/data/d-ibm3dx7008.txt", header = TRUE)
ibm = da$rtn
da[1:5,]
In [25]:
plot_acf(abs(ibm), lag.max = 100)
2-6 (Monthly Power Consumption)¶
- This is a multiplicative seasonal model.
- The seasonal prediction is very good.
- The data set in this problem bears a lot similarity with the JnJ quarterly earning data set in Sec. 2.8. After taking one regular difference, the seasonality becomes very obvious. This means that after taking one regular difference, we need to model the AR or MA pattern of the resulting time series.
- The logic of determining the orders in arima model from the various
pacfandacfplots can be as following.- The plots for
powshow strong serial correlation. So we need to take diff. - The plots for
diff(pow)show strong seasonality in theacf. So we need to take the seasonal diff. - The plots for
diff(pow, 12)show some autocorrelation in theacf. So we need to take both diff and seasonal diff. - The plots for
diff(diff(pow), 12)show some autocorrelation in theacf. So we need to apply MA model ondiff(diff(pow), 12). Thus this is a multiplicative seasonal model.
- The plots for
- Interesting observations.
- The
diff(log(pow))shows first strong momentum, second strong mean-reverting, and serial correlation inpacfandacf(b/c of seasonality). So the power-consumption time series doesn't satisfy "efficient market theory". So finanical instruments that are created around this time series have arbitrage opportunities. Seeplot_pacf_acfresults.
- The
In [502]:
da = read.table("../AFTS_sol/data/power6.txt", header = FALSE)
dim(da); da[1:5,]
In [507]:
par(bg = 'white')
pow = da[, 1]
pow_ts = ts(pow, frequency = 12, start = c(2000, 1))
# From the time plot, should immediately think that we can apply multiplicative seasonal model
plot(pow_ts, type = 'o', xlab = 'Year', ylab = 'Power', main = "Power Consumption")
plot(diff(pow_ts), type = 'o', xlab = 'Year', ylab = 'diff(Power)', main = "Power Consumption")
plot(diff(pow_ts, 12), type = 'o', xlab = 'Year', ylab = 'diff(Power, 12)', main = "Power Consumption")
plot(diff(diff(pow_ts), 12), type = 'o', xlab = 'Year', ylab = 'diff(diff(Power), 12)', main = "Power Consumption")
In [542]:
plot_pacf_acf(pow, freq = 12, lag.max = 48)
In [508]:
par(bg = 'white')
lg_pow = log(da[, 1])
lg_pow_ts = ts(lg_pow, frequency = 12, start = c(2000, 1))
# From the time plot, should immediately think that we can apply multiplicative seasonal model
plot(lg_pow_ts, type = 'o', xlab = 'Year', ylab = 'lg_power', main = "lg_power Consumption")
plot(diff(lg_pow_ts), type = 'o', xlab = 'Year', ylab = 'diff(lg_power)', main = "lg_power Consumption")
plot(diff(lg_pow_ts, 12), type = 'o', xlab = 'Year', ylab = 'diff(lg_power, 12)', main = "lg_power Consumption")
plot(diff(diff(lg_pow_ts), 12), type = 'o', xlab = 'Year', ylab = 'diff(diff(lg_power), 12)', main = "lg_power Consumption")
In [541]:
plot_pacf_acf(lg_pow, freq = 12, lag.max = 48)
In [ ]:
In [103]:
# The eacf is not helpful in determining the order of ARIMA model
pow_eacf_res <- perform_and_print_eacf(pow, ar.max = 25, ma.max = 12)
In [104]:
par(bg = 'white')
# FIXME
m1 = arima(pow, order = c(0, 1, 1), seasonal = list(order = c(0, 1, 1), period = 12))
m1
tsdiag(m1, gof = 36)
# predict(m1, 24)
In [111]:
pow_ts = ts(pow, frequency = 12, start = c(2010, 1))
npts = 24
eotr = length(pow_ts)-npts
h = npts
freq = 12
# The order and seasonal following in combination gives a multiplicative seasonal model.
order = c(0, 1, 1)
seasonal = list(order = c(0, 1, 1), period = 12)
fixed = NULL
pow_fc_res = plot_arima_forecast_fig(
da_ts=pow_ts, eotr=eotr, h=h, npts=npts, frequency=freq,
order=order, seasonal=seasonal, fixed=fixed, method='ML',
include.mean=T, transform.pars=TRUE,
main="Forecasts from ARIMA(0,1,1)(0,1,1)[12] for Power consumption",
xlab="Year", ylab="Power", ylim=c(9.6, 10.1)
)
In [44]:
pow_fc_tb = comb_forecast_res(pow_fc_res, pow_ts, eotr, freq)
pow_fc_tb
2-7 (Equal-Weighted Index)¶
In [512]:
require(xts)
require(sandwich)
require(lmtest)
require(forecast)
In [513]:
find("read.table")
Fit a time series model for CRSP equal-weighted index¶
- Data are daily simple returns.
- Test the effect of trading days using regression model without time series error.
- Add time series error model by using ARIMA model. The order of ARIMA can be seen from some acf and pacf plots; Also can look at eacf plots and tables. After getting a starting point for the
orderparameter ofarimafunction, we can trim the number of coefs by reducing element numbers inorder, until we have all coefs close to be statistically significant. - The logic of determining the order in
arimamodel is as following.- We see that the
m1$residualshave some autocorrelation inacf. We need to take diff. - After taking one diff, we see that the
diff(m1$residuals)have some autocorrelation in bothpacfandacf. We take seasonal diff. - After taking one seasonal diff, we see that
pacfhas some strong seasonality. Andacfhas some autocorrelation. Let's look at effect of both regular diff and seasonal diff. - After taking both regular and seasonal diff, we see that the autocorrelation in both
pacfandacfgetting stronger. We probably don't want to do this. - So, we need to apply seasonal diff and model the decaying autocorrelation in the
pacf. That means modeling the AR process.
- We see that the
- Interesting observations.
- The
log(ew_rtn)shows medium momentum and serial correlation inpacfat lag-1. Seeplot_pacf_acfresults.
- The
In [516]:
da = read.table("../AFTS_sol/data/d-ibm3dxwkdays8008.txt", header = TRUE)
da[1:5,]
In [517]:
ew = da$ew * 100
ew_ts = ts(ew, frequency = 252, start = c(1980, 1, 2))
plot(
xts(ew, order.by = as.Date(paste(da$year, da$mom, da$day, sep = '-'))),
type = 'l', main = '', xlab = 'date', ylab = 'ew'
)
In [539]:
plot_pacf_acf(ew, freq = 12, lag.max = 42)
In [519]:
lg_ew = log(1+da$ew)
lg_ew_ts = ts(lg_ew, frequency = 252, start = c(1980, 1, 2))
plot(
xts(lg_ew_ts, order.by = as.Date(paste(da$year, da$mom, da$day, sep = '-'))),
type = 'l', main = '', xlab = 'date', ylab = 'lg(1+ew)'
)
In [540]:
plot_pacf_acf(lg_ew, freq = 12, lag.max = 42)
- Trading day effects are all significant.
In [116]:
Mo = da$M
Tu = da$T
We = da$W
Th = da$R
m1 = lm(ew ~ Mo + Tu + We + Th)
summary(m1)
- Trading day effects are still significant. So HAC estimator doesn't change the conclusion of weekday effects.
In [17]:
find("NeweyWest")
In [117]:
coeftest(m1, NeweyWest(m1, lag = 12, prewhite = F))
- Residuals have serial correlation.
- Box-Ljung test tells us the null-hypothesis can be rejected.
In [118]:
Box.test(m1$residuals, lag = 12, type = 'Ljung')
- Residuals have some autocorrelations.
In [119]:
par(bg = 'white')
plot_time_fig(m1$residuals)
plot_time_fig(m1$residuals[1:200])
plot_time_fig(m1$residuals[(length(m1$residuals)-200):length(m1$residuals)])
In [122]:
par(mfrow = c(2, 1), bg = 'white')
pacf(m1$residuals)
acf(m1$residuals)
pacf(diff(m1$residuals))
acf(diff(m1$residuals))
pacf(diff(m1$residuals, 5))
acf(diff(m1$residuals, 5))
pacf(diff(diff(m1$residuals), 5))
acf(diff(diff(m1$residuals), 5))
- Fit an ARIMA model with time series error.
In [59]:
ew_eacf_res <- perform_and_print_eacf(ew, ar.max = 25, ma.max = 12)
In [56]:
par(bg = 'white')
m2 = arima(
ew, order = c(2, 0, 2),
seasonal = list(order = c(1, 0, 0), period = 5),
xreg = da[, 8:11]
)
m2
tsdiag(m2, gof = 20)
In [58]:
par(bg = 'white')
m3 = arima(
ew, order = c(2, 0, 3),
seasonal = list(order = c(1, 0, 0), period = 5),
xreg = da[, 8:11]
)
m3
tsdiag(m3, gof = 20)
In [24]:
par(bg = 'white')
m4 = arima(
ew, order = c(4, 0, 4),
seasonal = list(order = c(1, 0, 0), period = 5),
xreg = da[, 8:11]
)
m4
tsdiag(m4, gof = 20)
In [25]:
par(bg = 'white')
m4 = arima(
ew, order = c(3, 0, 3),
seasonal = list(order = c(1, 0, 0), period = 5),
xreg = da[, 8:11]
)
m4
tsdiag(m4, gof = 20)
In [126]:
par(bg = 'white')
m5 = arima(
ew, order = c(2, 0, 2),
# seasonal = list(order = c(0, 0, 1), period = 5),
xreg = da[, 8:11]
)
m5
tsdiag(m5, gof = 20)
In [125]:
par(bg = 'white')
m5 = arima(
ew, order = c(2, 0, 2),
seasonal = list(order = c(0, 0, 1), period = 5),
xreg = da[, 8:11]
)
m5
tsdiag(m5, gof = 20)
In [128]:
par(bg = 'white')
m5 = arima(
ew, order = c(2, 0, 2),
seasonal = list(order = c(1, 0, 0), period = 5),
xreg = da[, 8:11]
)
m5
tsdiag(m5, gof = 20)
In [127]:
par(bg = 'white')
m5 = arima(
ew, order = c(2, 0, 2),
seasonal = list(order = c(1, 0, 1), period = 5),
xreg = da[, 8:11]
)
m5
tsdiag(m5, gof = 20)
Forecast using the model¶
In [7]:
help(auto.arima)
Test auto.arima with a subset of time series data¶
In [17]:
da[6560:6570,]
In [19]:
sub_ew_ts = ts(ew[6565:length(ew)], frequency = 252, start = c(2006, 1, 3))
sub_da = da[6565:length(ew),]
length(sub_ew_ts)
In [23]:
sub_ts_fm <- auto.arima(
sub_ew_ts,
d = 0,
D = 0,
max.p = 2,
max.q = 2,
max.P = 1,
max.Q = 0,
max.order = 5,
# nmodels = 10,
seasonal = TRUE,
method = "ML",
allowmean = TRUE,
xreg = as.matrix(sub_da[, 8:11]),
# stepwise = FALSE,
# parallel = TRUE,
# num.cores = 6
)
sub_ts_fm
In [25]:
class(sub_ts_fm); attributes(sub_ts_fm)
In [29]:
h = 20
sub_ts_fc_res <- forecast(sub_ts_fm, h = h, xreg = as.matrix(sub_da[(dim(sub_da)[1]-h):(dim(sub_da)[1]), 8:11]))
summary(sub_ts_fc_res)
In [31]:
par(bg = 'white')
plot(sub_ts_fc_res)
In [32]:
sub_ts_fc_res$mean
In [33]:
sub_ts_fc_res$lower; sub_ts_fc_res$upper
Plot forecast results from auto.arima model with a subset of time series¶
In [41]:
help(stopifnot)
In [112]:
npts = 20
eotr = length(sub_ew_ts)-npts
h = npts
freq = 252
xreg = as.matrix(sub_da[, 8:11])
sub_ew_fc_res = plot_auto_arima_forecast_fig(
da_ts=sub_ew_ts, eotr=eotr, h=h, npts=npts, frequency=freq,
xreg=xreg,
main="Forecasts from ARIMA(1,1,0)\nfor CRSP Equal-Weighted Index",
xlab="Year", ylab="CRSP Equal-Weighted Index"# , ylim=c(5, 11)
)
In [49]:
sub_ew_fc_tb = comb_forecast_res(sub_ew_fc_res, sub_ew_ts, eotr, freq)
sub_ew_fc_tb
In [62]:
class(sub_ew_fc_res); attributes(sub_ew_fc_res);
class(sub_ew_fc_res$model); attributes(sub_ew_fc_res$model);
In [59]:
sub_ew_fc_res$model
All data takes too long time to run auto.arima¶
Dummy example for auto.arima¶
In [50]:
# Assuming 'forecast' package is installed and loaded
library(forecast)
# Simulate some time series data
set.seed(123)
ts_data <- ts(rnorm(100), frequency = 12, start = c(2000, 1))
# Simulate exogenous regressor data for model fitting
xreg_train <- matrix(rnorm(200), ncol = 2) # Two exogenous regressors
# Fit an ARIMA model with exogenous regressors
# fit <- auto.arima(ts_data, xreg = xreg_train)
fit <- auto.arima(ts_data, xreg = xreg_train)
# Simulate future values of the exogenous regressors for forecasting
xreg_future <- matrix(rnorm(24), ncol = 2) # Two exogenous regressors for 12 future periods
# Forecast using the model and future values of exogenous regressors
forecast_result <- forecast(fit, xreg = xreg_future)
# Plot the forecast
par(bg = 'white')
plot(forecast_result)
In [57]:
npts = 24
eotr = length(ts_data)-npts
h = npts
freq = 12
xreg = xreg_train
sim_data_fc_res = plot_auto_arima_forecast_fig(
da_ts=ts_data, eotr=eotr, h=h, npts=npts, frequency=freq,
xreg=xreg,
main="Forecasts from ARIMA(0,0,0) for Simulated Data",
xlab="Year", ylab="Simulated Data"# , ylim=c(5, 11)
)
2-8 (S&P)¶
- Data are daily simple returns.
- Interesting observations.
- The
diff(log(sp))shows medium mean-reverting and serial correlation inpacfat lag-2. Seeplot_pacf_acfresults.
- The
In [141]:
require(xts)
In [145]:
da = read.table("../AFTS_sol/data/d-ibm3dxwkdays8008.txt", header = TRUE)
dim(da); da[1:5,]
In [524]:
sp = da$sp * 100
sp_ts = ts(sp, frequency = 252, start = c(1980, 1, 2))
plot(xts(sp, order.by = as.Date(paste(da$year, da$mom, da$day, sep = '-'))),
type = 'l', main = '', xlab = 'date', ylab = 'sp')
lg_sp = log(1+da$sp)
lg_sp_ts = ts(lg_sp, frequency = 252, start = c(1980, 1, 2))
plot(xts(lg_sp, order.by = as.Date(paste(da$year, da$mom, da$day, sep = '-'))),
type = 'l', main = '', xlab = 'date', ylab = 'lg(sp)')
In [522]:
plot_pacf_acf(sp, freq = 12, lag.max = 42)
In [538]:
plot_pacf_acf(lg_sp, freq = 12, lag.max = 42)
2-8 (a)¶
- No Friday effect.
In [159]:
plot(
xts(sp, order.by = as.Date(paste(da$year, da$mom, da$day, sep = '-'))),
type = 'l', main = '', xlab = 'date', ylab = 'sp'
)
In [162]:
Mo = da$M
Tu = da$T
We = da$W
Th = da$R
Fr = da$F
sp_m1 = lm(sp ~ Mo + Tu + We + Th + Fr + 0)
summary(sp_m1);
In [163]:
Box.test(sp_m1$residuals, lag = 12, type = 'Ljung');
In [220]:
sp_m1 = lm(sp ~ 1 + T + W)
summary(sp_m1);
2-8 (b)¶
- There are some serial correlation in the residuals.
In [164]:
plot_pacf_acf(sp_m1$residuals, lag.max = 20, main = 'obs')
plot_pacf_acf(diff(sp_m1$residuals), lag.max = 20, main = 'diff(obs)')
plot_pacf_acf(diff(sp_m1$residuals, 5), lag.max = 20, main = "diff(obs, 5)")
plot_pacf_acf(diff(diff(sp_m1$residuals), 5), lag.max = 20, main = "diff(diff(obs), 5)")
In [169]:
par(bg = 'white')
sp_m2 = arima(
sp, order = c(0, 0, 1),
# seasonal = list(order = c(1, 0, 0), period = 5),
xreg = da[, 8:11]
)
sp_m2
tsdiag(sp_m2, gof = 20)
In [170]:
par(bg = 'white')
sp_m2 = arima(
sp, order = c(1, 0, 0),
# seasonal = list(order = c(1, 0, 0), period = 5),
xreg = da[, 8:11]
)
sp_m2
tsdiag(sp_m2, gof = 20)
In [171]:
par(bg = 'white')
sp_m2 = arima(
sp, order = c(2, 0, 0),
# seasonal = list(order = c(1, 0, 0), period = 5),
xreg = da[, 8:11]
)
sp_m2
tsdiag(sp_m2, gof = 20)
In [172]:
par(bg = 'white')
sp_m2 = arima(
sp, order = c(2, 0, 0),
seasonal = list(order = c(1, 0, 0), period = 5),
xreg = da[, 8:11]
)
sp_m2
tsdiag(sp_m2, gof = 20)
In [174]:
par(bg = 'white')
sp_m2 = arima(
sp, order = c(3, 0, 1),
xreg = da[, 8:11]
)
sp_m2
tsdiag(sp_m2, gof = 20)
In [195]:
par(bg = 'white')
sp_m2 = arima(
sp, order = c(4, 0, 0),
xreg = da[, 8:11]
)
sp_m2
tsdiag(sp_m2, gof = 20)
In [196]:
par(bg = 'white')
sp_m2 = arima(
sp, order = c(4, 0, 0),
fixed = c(NA, NA, 0, NA, 0, 0, 0, 0, 0),
xreg = da[, 8:11]
)
sp_m2
tsdiag(sp_m2, gof = 20)
In [198]:
par(bg = 'white')
sp_m2 = arima(
sp, order = c(5, 0, 0),
fixed = c(NA, NA, 0, NA, 0, 0, 0, 0, 0, 0),
xreg = da[, 8:11]
)
sp_m2
tsdiag(sp_m2, gof = 20)
In [ ]:
In [200]:
par(bg = 'white')
sp_m3 = arima(
sp, order = c(2, 0, 0),
xreg = da[, 8:11]
)
sp_m3
tsdiag(sp_m3, gof = 20)
In [201]:
par(bg = 'white')
sp_m3 = arima(
sp, order = c(2, 0, 2),
xreg = da[, 8:11]
)
sp_m3
tsdiag(sp_m3, gof = 20)
In [203]:
par(bg = 'white')
sp_m3 = arima(
sp, order = c(2, 0, 2),
fixed = c(0, NA, 0, NA, 0, 0, 0, 0, 0),
xreg = da[, 8:11]
)
sp_m3
tsdiag(sp_m3, gof = 20)
In [ ]:
In [ ]:
In [213]:
sp_m4 = arima(sp, order = c(2, 0, 2), xreg = da[, 8:12], include.mean = FALSE)
sp_m4;
par(bg = 'white')
tsdiag(sp_m4, gof = 20)
- Adding T and W, b/c they show significance in linear model.
In [215]:
sp_m4 = arima(sp, order = c(2, 0, 2), xreg = da[, c(9, 10, 12)], include.mean = FALSE,
fixed = c(0, NA, 0, NA, NA, NA, NA))
sp_m4;
par(bg = 'white')
tsdiag(sp_m4, gof = 20)
Box.test(sp_m4$residuals, lag = 10, type = 'Ljung')
In [ ]:
2-9 (IBM)¶
- Data are daily simple returns.
- There is almost no serial correlation. See
plot_pacf_acfresults.
In [217]:
require(xts)
In [525]:
da = read.table("../AFTS_sol/data/d-ibm3dxwkdays8008.txt", header = TRUE)
da[1:5,]
In [526]:
ibm = da$ibm * 100
plot(xts(ibm, order.by = as.Date(paste(da$year, da$mom, da$day, sep = '-'))),
type = 'l', main = '', xlab = 'date', ylab = 'ibm')
lg_ibm = log(1+da$ibm)
plot(xts(lg_ibm, order.by = as.Date(paste(da$year, da$mom, da$day, sep = '-'))),
type = 'l', main = '', xlab = 'date', ylab = 'lg(ibm)')
In [535]:
plot_pacf_acf(ibm, freq = 12, lag.max = 42)
In [536]:
plot_pacf_acf(lg_ibm, freq = 12, lag.max = 42)
2-9 (a)¶
- No Friday effect.
In [223]:
Mo = da$M
Tu = da$T
We = da$W
Th = da$R
Fr = da$F
ibm_m1 = lm(ibm ~ Mo + Tu + We + Th + Fr + 0)
summary(ibm_m1)
In [225]:
Box.test(ibm_m1$residuals, lag = 12, type = 'Ljung')
2-9 (b)¶
- There is no much serial correlation in the residuals.
In [231]:
lag.max = 20
par(mfrow = c(2, 1), bg = 'white')
pacf(ibm_m1$residuals, lag.max =lag.max)
acf(ibm_m1$residuals, lag.max =lag.max)
pacf(diff(ibm_m1$residuals), lag.max =lag.max)
acf(diff(ibm_m1$residuals), lag.max =lag.max)
pacf(diff(ibm_m1$residuals, 5), lag.max =lag.max)
acf(diff(ibm_m1$residuals, 5), lag.max =lag.max)
pacf(diff(diff(ibm_m1$residuals), 5), lag.max =lag.max)
acf(diff(diff(ibm_m1$residuals), 5), lag.max =lag.max)
2-9 (c)¶
In [228]:
ibm_m2 = arima(ibm, order = c(0, 0, 1), xreg = da[, 8:12], include.mean = FALSE)
ibm_m2;
par(bg = 'white')
tsdiag(ibm_m2, gof = 35)
In [229]:
Box.test(ibm_m2$residuals, lag = 12, type = 'Ljung')
In [234]:
ibm_m2 = arima(ibm, order = c(0, 0, 1), xreg = da[, c(8, 9)], include.mean = TRUE)
ibm_m2;
par(bg = 'white')
tsdiag(ibm_m2, gof = 35)
In [233]:
ibm_m2 = arima(ibm, order = c(0, 0, 1), xreg = da[, c(8, 9)], include.mean = FALSE)
ibm_m2;
par(bg = 'white')
tsdiag(ibm_m2, gof = 35)
2-10 (Aaa and Baa Bond Yields)¶
- pp10: Test statistics for skewness and kurtosis.
- They are all significant.
In [237]:
require(fBasics)
In [235]:
da1 = read.table("../AFTS_sol/data/w-Aaa.txt", header = FALSE)
da2 = read.table("../AFTS_sol/data/w-Baa.txt", header = FALSE)
da1[1:5,]; da2[1:5,];
In [238]:
Aaa = da1[, 4]
Baa = da2[, 4]
apply(cbind(Aaa, Baa), 2, basicStats)
ts = skewness(Aaa) / sqrt(6 / length(Aaa))
cat("ts =", ts, "\n")
ps = (1 - pnorm(ts)) * 2
cat("ps =", ps, "\n")
tk = kurtosis(Aaa) / sqrt(24 / length(Aaa))
cat("tk =", tk, "\n")
pk = (1 - pnorm(tk)) * 2
cat("pk =", pk, "\n")
ts = skewness(Baa) / sqrt(6 / length(Baa))
cat("ts =", ts, "\n")
ps = (1 - pnorm(ts)) * 2
cat("ps =", ps, "\n")
tk = kurtosis(Baa) / sqrt(24 / length(Baa))
cat("tk =", tk, "\n")
pk = (1 - pnorm(tk)) * 2
cat("pk =", pk, "\n")
2-11 (Aaa and Baa Bond Yields)¶
- Data are weekly Bond Yields.
- Interesting observation.
- The
diff(aaa)ordiff(log(aaa))shows medium to strong momentum at lag-1 inpacfandacf, but not in longer lags. Seeplot_pacf_acfresults.
- The
In [345]:
da1 = read.table("../AFTS_sol/data/w-Aaa.txt", header = FALSE)
da1[1:5,];
In [346]:
aaa = da1[,4]
aaa_ts = ts(aaa, frequency = 52, start = c(1962, 1))
lg_aaa = log(aaa)
lg_aaa_ts = ts(lg_aaa, frequency = 52, start = c(1962, 1))
Fitting solution¶
In [253]:
sub_sz = 100
plot_time_fig(lg_aaa_ts, main = "AAA Bond", xlab = "Year", ylab = "lg(Bond Weekly Yield)")
plot_time_fig(lg_aaa_ts[1:sub_sz], main = "AAA Bond", xlab = "Year", ylab = "lg(Bond Weekly Yield)")
plot_time_fig(lg_aaa_ts[(length(lg_aaa_ts)-sub_sz):length(lg_aaa_ts)], main = "AAA Bond", xlab = "Year", ylab = "lg(Bond Weekly Yield)")
In [243]:
sub_sz = 100
plot_time_fig(aaa_ts, main = "AAA Bond", xlab = "Year", ylab = "Bond Weekly Yield")
plot_time_fig(aaa_ts[1:sub_sz], main = "AAA Bond", xlab = "Year", ylab = "Bond Weekly Yield")
plot_time_fig(aaa_ts[(length(aaa_ts)-sub_sz):length(aaa_ts)], main = "AAA Bond", xlab = "Year", ylab = "Bond Weekly Yield")
In [250]:
# lag.max = 52
# par(mfrow = c(2, 1), bg = 'white')
# pacf(aaa, lag.max = lag.max)
# acf(aaa, lag.max = lag.max)
# pacf(diff(aaa), lag.max = lag.max)
# acf(diff(aaa), lag.max = lag.max)
# pacf(diff(aaa, 4), lag.max = lag.max)
# acf(diff(aaa, 4), lag.max = lag.max)
# pacf(diff(aaa, 13), lag.max = lag.max)
# acf(diff(aaa, 13), lag.max = lag.max)
# pacf(diff(aaa, 26), lag.max = lag.max)
# acf(diff(aaa, 26), lag.max = lag.max)
# pacf(diff(diff(aaa), 26), lag.max = lag.max)
# acf(diff(diff(aaa), 26), lag.max = lag.max)
In [543]:
plot_pacf_acf(aaa, freq = 26, lag.max = 52)
In [544]:
plot_pacf_acf(lg_aaa, freq = 26, lag.max = 52)
In [ ]:
In [256]:
par(bg = 'white')
aaa_m1 = arima(
aaa, order = c(0, 1, 1),
seasonal = list(order = c(0, 1, 1), period = 26)
)
aaa_m1
tsdiag(aaa_m1, gof = 52)
In [259]:
par(bg = 'white')
aaa_m2 = arima(
aaa, order = c(1, 1, 0),
seasonal = list(order = c(1, 1, 0), period = 26)
)
aaa_m2
tsdiag(aaa_m2, gof = 52)
In [258]:
par(bg = 'white')
aaa_m2 = arima(
aaa, order = c(1, 1, 1),
seasonal = list(order = c(1, 1, 1), period = 26)
)
aaa_m2
tsdiag(aaa_m2, gof = 52)
In [260]:
par(bg = 'white')
aaa_m2 = arima(
aaa, order = c(1, 1, 1),
seasonal = list(order = c(0, 1, 1), period = 26)
)
aaa_m2
tsdiag(aaa_m2, gof = 52)
In [261]:
par(bg = 'white')
aaa_m2 = arima(
aaa, order = c(2, 1, 1),
seasonal = list(order = c(0, 1, 1), period = 26)
)
aaa_m2
tsdiag(aaa_m2, gof = 52)
In [262]:
par(bg = 'white')
aaa_m2 = arima(
aaa, order = c(3, 1, 1),
seasonal = list(order = c(0, 1, 1), period = 26)
)
aaa_m2
tsdiag(aaa_m2, gof = 52)
In [263]:
par(bg = 'white')
aaa_m2 = arima(
aaa, order = c(4, 1, 1),
seasonal = list(order = c(0, 1, 1), period = 26)
)
aaa_m2
tsdiag(aaa_m2, gof = 52)
In [265]:
par(bg = 'white')
aaa_m2 = arima(
aaa, order = c(9, 1, 1),
seasonal = list(order = c(0, 1, 1), period = 26)
)
aaa_m2
tsdiag(aaa_m2, gof = 52)
In [266]:
par(bg = 'white')
aaa_m2 = arima(
aaa, order = c(9, 1, 1),
seasonal = list(order = c(0, 1, 1), period = 26),
fixed = c(NA, 0, NA, 0, NA, NA, NA, NA, NA, 0, NA)
)
aaa_m2
tsdiag(aaa_m2, gof = 52)
In [267]:
par(bg = 'white')
aaa_m2 = arima(
aaa, order = c(9, 1, 1),
seasonal = list(order = c(0, 1, 1), period = 26),
fixed = c(NA, 0, NA, 0, NA, 0, NA, 0, NA, 0, NA)
)
aaa_m2
tsdiag(aaa_m2, gof = 52)
In [273]:
par(bg = 'white')
aaa_m2 = arima(aaa, order = c(9, 1, 1))
aaa_m2
tsdiag(aaa_m2, gof = 52)
In [276]:
par(bg = 'white')
aaa_m2 = arima(
aaa, order = c(9, 1, 0),
fixed = c(NA, 0, NA, 0, 0, 0, NA, 0, NA)
)
aaa_m2
tsdiag(aaa_m2, gof = 52)
In [ ]:
Alternative solution¶
In [268]:
plot(
xts(aaa, order.by = as.Date(paste(da1[, 1], da1[, 2], da1[, 3], sep = '-'))),
type = 'l', main = '', xlab = 'date', ylab = 'Aaa'
)
In [270]:
cat("order =", ar(diff(aaa), method = 'mle')$order, "\n")
In [271]:
aaa_al_m1 = arima(aaa, order = c(9, 1, 0))
aaa_al_m1;
par(bg = 'white')
tsdiag(aaa_al_m1, gof = 52)
In [291]:
aaa_al_m1 = arima(
aaa, order = c(9, 1, 0),
seasonal = list(order = c(1, 0, 0), period = 13)
)
aaa_al_m1;
par(bg = 'white')
tsdiag(aaa_al_m1, gof = 52)
In [277]:
aaa_al_m2 = arima(aaa, order = c(2, 1, 3))
aaa_al_m2;
par(bg = 'white')
tsdiag(aaa_al_m2, gof = 52)
In [288]:
aaa_al_m2 = arima(
aaa, order = c(2, 1, 3),
seasonal = list(order = c(1, 0, 0), period = 26)
)
aaa_al_m2;
par(bg = 'white')
tsdiag(aaa_al_m2, gof = 52)
Forecast¶
In [405]:
npts = 26
eotr = length(aaa_ts)-npts
h = npts
freq = 52
order = c(9,1,0)
fixed = NULL
seasonal = NULL
aaa_fc_res = plot_arima_forecast_fig(
da_ts=aaa_ts, eotr=eotr, h=h, npts=npts, frequency=freq,
order=order, seasonal=seasonal, fixed=fixed, method='ML',
include.mean=TRUE, transform.pars=TRUE,
main="Forecasts from ARIMA(9,1,0) for Aaa weekly bond yields",
xlab="Year", ylab="Aaa weekly bond yields", ylim=c(4, 8)
)
In [283]:
# aaa_fc_tb = comb_forecast_res(aaa_fc_res, aaa_ts, eotr, freq)
# aaa_fc_tb
In [495]:
aaa_fc_tb = comb_forecast_res(aaa_fc_res, aaa_ts, eotr, freq)
aaa_fc_res$model['coef'];
sqrt(diag(aaa_fc_res$model[['var.coef']]));
aaa_fc_tb
In [402]:
npts = 26
eotr = length(aaa_ts)-npts
h = npts
freq = 52
order = c(2,1,3)
fixed = NULL
seasonal = NULL
aaa_fc_res_1 = plot_arima_forecast_fig(
da_ts=aaa_ts, eotr=eotr, h=h, npts=npts, frequency=freq,
order=order, seasonal=seasonal, fixed=fixed, method='ML',
include.mean=TRUE, transform.pars=TRUE,
main="Forecasts from ARIMA(2,1,3) for Aaa weekly bond yields",
xlab="Year", ylab="Aaa weekly bond yields", ylim=c(4, 8)
)
In [285]:
# aaa_fc_tb_1 = comb_forecast_res(aaa_fc_res_1, aaa_ts, eotr, freq)
# aaa_fc_tb_1
In [498]:
aaa_fc_tb_1 = comb_forecast_res(aaa_fc_res_1, aaa_ts, eotr, freq)
aaa_fc_res_1$model['coef'];
sqrt(diag(aaa_fc_res_1$model[['var.coef']]));
aaa_fc_tb_1
In [373]:
help(arima)
In [499]:
npts = 26
eotr = length(aaa_ts)-npts
h = npts
freq = 52
order = c(2,1,3)
fixed = NULL
seasonal = list(order = c(0, 0, 1), period = 13)
aaa_fc_res_2 = plot_arima_forecast_fig(
da_ts=aaa_ts, eotr=eotr, h=h, npts=npts, frequency=freq,
order=order, seasonal=seasonal, fixed=fixed, method='ML',
include.mean=TRUE, transform.pars=TRUE,
main=NULL, #"Forecasts from ARIMA(2,1,3)(0,0,1)[13] for Aaa weekly bond yields",
xlab="Year", ylab="Aaa weekly bond yields", ylim=c(4, 8)
)
In [500]:
aaa_fc_tb_2 = comb_forecast_res(aaa_fc_res_2, aaa_ts, eotr, freq)
aaa_fc_res_2$model['coef'];
sqrt(diag(aaa_fc_res_2$model[['var.coef']]));
aaa_fc_tb_2
In [ ]:
2-12 (Aaa and Baa Bond Yields)¶
In [299]:
da1 = read.table("../AFTS_sol/data/w-Aaa.txt", header = FALSE)
da2 = read.table("../AFTS_sol/data/w-Baa.txt", header = FALSE)
da1[1:5,]; da2[1:5,];
In [302]:
freq = 52
start = c(1962, 1)
aaa = da1[,4]
aaa_ts = ts(aaa, frequency = freq, start = start)
bbb = da2[,4]
bbb_ts = ts(bbb, frequency = freq, start = start)
aaac = diff(aaa)
bbbc = diff(bbb)
aaac_ts = ts(aaac, frequency = freq, start = start)
bbbc_ts = ts(bbbc, frequency = freq, start = start)
In [316]:
apply(cbind(aaac, bbbc), 2, basicStats)
Exploration¶
In [301]:
par(bg = 'white')
plot(aaa_ts, main = "Weekly Aaa or Bbb Bond Yield", xlab = "Year", ylab = "Yield")
lines(bbb_ts, lty = 2)
legend(
"topleft",
legend = c("Aaa", "Bbb"),
lty = c(1, 2)
)
In [304]:
par(bg = 'white')
plot(aaac_ts, main = "Weekly Aaa Bond Yield", xlab = "Year", ylab = "Yield")
plot(bbbc_ts, main = "Weekly Bbb Bond Yield", xlab = "Year", ylab = "Yield")
In [305]:
par(bg = "white")
plot(x = aaa, y = bbb, main = "Aaa vs Bbb", xlab = "Aaa Yield", ylab = "Bbb Yield")
plot(x = aaac, y = bbbc, main = "Aaa Change vs Bbb Change", xlab = "diff(Aaa Yield)", ylab = "diff(Bbb Yield)")
Fitting¶
- Strong serial correlation in ACF of residuals, indicating that there is a unit-root process.
- Also conducted ADF unit-root tests. The results show that null-hypothesis cannot be rejected.
- Null-hypothesis is that the time series has unit-root.
- Box-Ljung test shows that numm-hypothesis can be rejected.
- Null-hypothesis is that the time series at certain lag has no autocorrelation.
- After taking one diff, the autocorrelation is largely gone. But Box-Ljung test is still significant. We can model the residuals as a time series.
Fitting the bond yields¶
In [306]:
bond_m1 = lm(aaa~bbb)
summary(bond_m1)
In [308]:
par(bg = 'white')
resi_ts = ts(bond_m1$residuals, frequency = freq, start = start)
plot(resi_ts, type = 'l', main = "Simple Linear Model", xlab = "Year")
lag.max = 36
par(mfrow = c(2,1))
pacf(bond_m1$residuals, lag.max = lag.max, main = "Residuals PACF Plot")
acf(bond_m1$residuals, lag.max = lag.max, main = "Residuals ACF Plot")
In [309]:
adf_test_res_a = adfTest(aaa, lags = 12, type = c("ct"))
adf_test_res_b = adfTest(bbb, lags = 12, type = c("ct")) # 2.7.3. Trend-Stationary Time-Series
adf_test_res_a; adf_test_res_b
In [310]:
Box.test(bond_m1$residuals, lag = 12, type = 'Ljung')
Fitting the bond yield changes¶
In [311]:
bond_m2 = lm(aaac~-1+bbbc) # No intercept
summary(bond_m2)
In [312]:
par(bg = 'white')
resi_ts = ts(bond_m2$residuals, frequency = freq, start = start)
plot(resi_ts, type = 'l', main = "Simple Linear Model", xlab = "Year")
lag.max = 36
par(mfrow = c(2,1))
pacf(bond_m2$residuals, lag.max = lag.max, main = "Residuals PACF Plot")
acf(bond_m2$residuals, lag.max = lag.max, main = "Residuals ACF Plot")
In [314]:
adf_test_res_ac = adfTest(aaac, lags = 12, type = c("ct"))
adf_test_res_bc = adfTest(bbbc, lags = 12, type = c("ct")) # 2.7.3. Trend-Stationary Time-Series
adf_test_res_ac; adf_test_res_bc
In [313]:
Box.test(bond_m2$residuals, lag = 12, type = 'Ljung')
Fitting a time series residuals¶
In [587]:
bond_m3 = arima(aaac, order = c(0,0,1), xreg = bbbc, include.mean = FALSE)
# rsq = 1 - sum(m3$residuals^2)/sum(c3^2-mean(c3)) # B/c in the model we set `include.mean = F`, so here when evaluating R-squared, we don't include the mean either.
rsq = 1 - sum(bond_m3$residuals^2)/sum(aaac^2)
summary(bond_m3); rsq;
par(bg = 'white')
tsdiag(bond_m3, gof = 36)
In [588]:
bond_m3 = arima(aaac, order = c(1,0,0), xreg = bbbc, include.mean = FALSE)
# rsq = 1 - sum(m3$residuals^2)/sum(c3^2-mean(c3)) # B/c in the model we set `include.mean = F`, so here when evaluating R-squared, we don't include the mean either.
rsq = 1 - sum(bond_m3$residuals^2)/sum(aaac^2)
summary(bond_m3); rsq;
par(bg = 'white')
tsdiag(bond_m3, gof = 36)
In [589]:
bond_m3 = arima(aaac, order = c(1,0,1), xreg = bbbc, include.mean = FALSE)
# rsq = 1 - sum(m3$residuals^2)/sum(c3^2-mean(c3)) # B/c in the model we set `include.mean = F`, so here when evaluating R-squared, we don't include the mean either.
rsq = 1 - sum(bond_m3$residuals^2)/sum(aaac^2)
summary(bond_m3); rsq;
par(bg = 'white')
tsdiag(bond_m3, gof = 36)
In [ ]:
In [324]:
cat("order =", ar(bond_m2$residuals, method = 'mle')$order, "\n")
In [579]:
bond_m3 = arima(aaac, order = c(12,0,0), xreg = bbbc, include.mean = FALSE)
# rsq = 1 - sum(m3$residuals^2)/sum(c3^2-mean(c3)) # B/c in the model we set `include.mean = F`, so here when evaluating R-squared, we don't include the mean either.
rsq = 1 - sum(bond_m3$residuals^2)/sum(aaac^2)
summary(bond_m3); rsq;
par(bg = 'white')
tsdiag(bond_m3, gof = 36)
In [580]:
bond_m3 = arima(aaac, order = c(2,0,0), xreg = bbbc, include.mean = FALSE)
# rsq = 1 - sum(m3$residuals^2)/sum(c3^2-mean(c3)) # B/c in the model we set `include.mean = F`, so here when evaluating R-squared, we don't include the mean either.
rsq = 1 - sum(bond_m3$residuals^2)/sum(aaac^2)
summary(bond_m3); rsq;
par(bg = 'white')
tsdiag(bond_m3, gof = 36)
In [581]:
bond_m3 = arima(aaac, order = c(3,0,0), xreg = bbbc, include.mean = FALSE)
# rsq = 1 - sum(m3$residuals^2)/sum(c3^2-mean(c3)) # B/c in the model we set `include.mean = F`, so here when evaluating R-squared, we don't include the mean either.
rsq = 1 - sum(bond_m3$residuals^2)/sum(aaac^2)
summary(bond_m3); rsq;
par(bg = 'white')
tsdiag(bond_m3, gof = 36)
In [582]:
bond_m3 = arima(aaac, order = c(4,0,0), xreg = bbbc, include.mean = FALSE)
# rsq = 1 - sum(m3$residuals^2)/sum(c3^2-mean(c3)) # B/c in the model we set `include.mean = F`, so here when evaluating R-squared, we don't include the mean either.
rsq = 1 - sum(bond_m3$residuals^2)/sum(aaac^2)
summary(bond_m3); rsq;
par(bg = 'white')
tsdiag(bond_m3, gof = 36)
In [583]:
bond_m3 = arima(aaac, order = c(5,0,0), xreg = bbbc, include.mean = FALSE)
# rsq = 1 - sum(m3$residuals^2)/sum(c3^2-mean(c3)) # B/c in the model we set `include.mean = F`, so here when evaluating R-squared, we don't include the mean either.
rsq = 1 - sum(bond_m3$residuals^2)/sum(aaac^2)
summary(bond_m3); rsq;
par(bg = 'white')
tsdiag(bond_m3, gof = 36)
In [584]:
bond_m3 = arima(aaac, order = c(8,0,0), xreg = bbbc, include.mean = FALSE)
# rsq = 1 - sum(m3$residuals^2)/sum(c3^2-mean(c3)) # B/c in the model we set `include.mean = F`, so here when evaluating R-squared, we don't include the mean either.
rsq = 1 - sum(bond_m3$residuals^2)/sum(aaac^2)
summary(bond_m3); rsq;
par(bg = 'white')
tsdiag(bond_m3, gof = 36)
In [585]:
bond_m3 = arima(aaac, order = c(2,0,2), xreg = bbbc, include.mean = FALSE)
# rsq = 1 - sum(m3$residuals^2)/sum(c3^2-mean(c3)) # B/c in the model we set `include.mean = F`, so here when evaluating R-squared, we don't include the mean either.
rsq = 1 - sum(bond_m3$residuals^2)/sum(aaac^2)
summary(bond_m3); rsq;
par(bg = 'white')
tsdiag(bond_m3, gof = 36)
In [586]:
bond_m3 = arima(
aaac, order = c(2,0,2),
fixed = c(0, NA, NA, NA, NA),
xreg = bbbc, include.mean = FALSE
)
# rsq = 1 - sum(m3$residuals^2)/sum(c3^2-mean(c3)) # B/c in the model we set `include.mean = F`, so here when evaluating R-squared, we don't include the mean either.
rsq = 1 - sum(bond_m3$residuals^2)/sum(aaac^2)
summary(bond_m3); rsq;
par(bg = 'white')
tsdiag(bond_m3, gof = 36)
Forecast¶
In [332]:
npts = 26
eotr = length(aaac_ts)-npts
h = npts
freq = 52
xreg = as.matrix(bbbc)
aaac_fc_res = plot_auto_arima_forecast_fig(
da_ts=aaac_ts, eotr=eotr, h=h, npts=npts, frequency=freq,
xreg=xreg,
main=NULL, # "Forecasts from ARIMA(2,0,2)\nfor CRSP Equal-Weighted Index",
xlab="Year", ylab="Aaa Bond Yield (using Baa Bond Yield)", # , ylim=c(5, 11)
d = 0,
D = 0,
max.p = 3,
max.q = 4,
max.P = 0,
max.Q = 0,
max.order = 5,
seasonal = FALSE,
method = "ML",
allowmean = FALSE,
stepwise = FALSE,
parallel = TRUE,
num.cores = 2
)
In [333]:
aaac_fc_tb = comb_forecast_res(aaac_fc_res, aaac_ts, eotr, freq)
aaac_fc_tb
In [342]:
par(bg = 'white')
tsdiag(aaac_fc_res$model, gof.lag = 36)
In [ ]:
2-13 (CRSP Equal-Weighted Index)¶
- The data are monthly log-return.
- Interesting observations.
- The
ew_lrtnhas medium momentum at lag-1, and no autocorrelation for longer lags.
- The
In [453]:
library(lmtest)
In [454]:
da = read.table("../AFTS_sol/data/m-ew6299.txt", header = FALSE)
ew_lrtn = da[,1]
ew_lrtn_ts = ts(ew_lrtn, frequency = 12, start = c(1962, 1))
ew_prr = cumsum(ew_lrtn)
ew_prr_ts = ts(ew_prr, frequency = 12, start = c(1962, 1))
dim(da); da[1:5, ]
In [455]:
start(ew_prr_ts)
In [659]:
t_offset = c(10, 0)
plot_time_fig(ew_prr_ts, main = "CRSP Equal-Weighted Index", xlab = "Year", ylab = "lg(pr_ratio)")
plot_time_fig(
window(ew_prr_ts, start = start(ew_prr_ts), end = start(ew_prr_ts) + t_offset),
main = "CRSP Equal-Weighted Index", xlab = "Year", ylab = "lg(pr_ratio)"
)
plot_time_fig(
window(ew_prr_ts, start = end(ew_prr_ts) - t_offset, end = end(ew_prr_ts)),
main = "CRSP Equal-Weighted Index", xlab = "Year", ylab = "lg(pr_ratio)"
)
In [660]:
t_offset = c(10, 0)
plot_time_fig(ew_lrtn_ts, main = "CRSP Equal-Weighted Index", xlab = "Year", ylab = "lg(rtn)")
plot_time_fig(
window(ew_lrtn_ts, start = start(ew_lrtn_ts), end = start(ew_lrtn_ts) + t_offset),
main = "CRSP Equal-Weighted Index", xlab = "Year", ylab = "lg(pr_ratio)"
)
plot_time_fig(
window(ew_lrtn_ts, start = end(ew_lrtn_ts) - t_offset, end = end(ew_lrtn_ts)),
main = "CRSP Equal-Weighted Index", xlab = "Year", ylab = "lg(pr_ratio)"
)
In [573]:
plot_pacf_acf(ew_lrtn, freq = 12, lag.max = 36)
2-13 (a)¶
- The
AR(1, 0, 0)looks adequate.
In [577]:
cat("order =", ar(ew_lrtn, method = 'mle')$order, "\n")
In [639]:
ew_lrtn_m1 = arima(ew_lrtn_ts, order = c(1, 0, 0), method = 'ML')
rsq = 1 - sum(ew_lrtn_m1$residuals^2)/sum(ew_lrtn^2)
summary(ew_lrtn_m1); rsq;
par(bg = 'white')
tsdiag1(ew_lrtn_m1, gof.lag = 36)
In [646]:
ew_lrtn_m2 = arima(ew_lrtn_ts, order = c(2, 0, 0), method = 'ML')
rsq = 1 - sum(ew_lrtn_m2$residuals^2)/sum(ew_lrtn^2)
summary(ew_lrtn_m2); rsq;
par(bg = 'white')
tsdiag1(ew_lrtn_m2, gof = 36)
2-13 (b)¶
- The
MA(0, 0, 1)looks adequate.
In [737]:
ew_lrtn_m3 = arima(ew_lrtn_ts, order = c(0, 0, 1), method = 'ML')
rsq = 1 - sum(ew_lrtn_m3$residuals^2)/sum(ew_lrtn^2)
summary(ew_lrtn_m3); rsq;
par(bg = 'white')
tsdiag1(ew_lrtn_m3, gof = 36)
# likelihood-ratio test
llike_rtest_pv <- pchisq(
q = 2*(ew_lrtn_m3$loglik-ew_lrtn_m1$loglik),
df = length(ew_lrtn_m3$coef)-length(ew_lrtn_m1$coef),
lower.tail = FALSE
)
# cat("Likelihood-Ratio test p-value:", llike_rtest_pv, "")
lrt_res <- lrtest(ew_lrtn_m1, ew_lrtn_m3)
sprintf(
"Likelihood-Ratio test p-value: %f (lrtest: %d; %f; %f)",
llike_rtest_pv,
lrt_res$Df[2],
lrt_res$Chisq[2],
lrt_res$`Pr(>Chisq)`[2]
)
In [648]:
ew_lrtn_m4 = arima(ew_lrtn_ts, order = c(0, 0, 2), method = 'ML')
rsq = 1 - sum(ew_lrtn_m4$residuals^2)/sum(ew_lrtn^2)
summary(ew_lrtn_m4); rsq;
par(bg = 'white')
tsdiag1(ew_lrtn_m4, gof = 36)
In [650]:
ew_lrtn_m5 = arima(ew_lrtn_ts, order = c(1, 0, 1), method = 'ML')
rsq = 1 - sum(ew_lrtn_m5$residuals^2)/sum(ew_lrtn^2)
summary(ew_lrtn_m5); rsq;
par(bg = 'white')
tsdiag1(ew_lrtn_m5, gof = 36)
In [736]:
ew_lrtn_m6 = arima(ew_lrtn_ts, order = c(1, 0, 0), seasonal = c(1, 0, 0), method = 'ML')
rsq = 1 - sum(ew_lrtn_m6$residuals^2)/sum(ew_lrtn^2)
summary(ew_lrtn_m6); rsq;
par(bg = 'white')
tsdiag1(ew_lrtn_m6, gof = 36)
# likelihood-ratio test
llike_rtest_pv <- pchisq(
q = 2*(ew_lrtn_m6$loglik-ew_lrtn_m1$loglik),
df = length(ew_lrtn_m6$coef)-length(ew_lrtn_m1$coef),
lower.tail = FALSE
)
# cat("Likelihood-Ratio test p-value:", llike_rtest_pv, "")
lrt_res <- lrtest(ew_lrtn_m1, ew_lrtn_m6)
sprintf(
"Likelihood-Ratio test p-value: %f (lrtest: %d; %f; %f)",
llike_rtest_pv,
lrt_res$Df[2],
lrt_res$Chisq[2],
lrt_res$`Pr(>Chisq)`[2]
)
In [735]:
ew_lrtn_m7 = arima(ew_lrtn_ts, order = c(0, 0, 1), seasonal = c(0, 0, 1), method = 'ML')
rsq = 1 - sum(ew_lrtn_m7$residuals^2)/sum(ew_lrtn^2)
summary(ew_lrtn_m7); rsq;
par(bg = 'white')
tsdiag1(ew_lrtn_m7, gof = 36)
# likelihood-ratio test
llike_rtest_pv <- pchisq(
q = 2*(ew_lrtn_m7$loglik-ew_lrtn_m3$loglik),
df = length(ew_lrtn_m7$coef)-length(ew_lrtn_m3$coef),
lower.tail = FALSE
)
# cat("Likelihood-Ratio test p-value:", llike_rtest_pv, "")
lrt_res <- lrtest(ew_lrtn_m3, ew_lrtn_m7)
sprintf(
"Likelihood-Ratio test p-value: %f (lrtest: %d; %f; %f)",
llike_rtest_pv,
lrt_res$Df[2],
lrt_res$Chisq[2],
lrt_res$`Pr(>Chisq)`[2]
)
In [744]:
ew_lrtn_m8 = arima(ew_lrtn_ts, order = c(1, 0, 0), seasonal = c(1, 0, 1), method = 'ML')
rsq = 1 - sum(ew_lrtn_m8$residuals^2)/sum(ew_lrtn^2)
summary(ew_lrtn_m8); rsq;
par(bg = 'white')
tsdiag1(ew_lrtn_m8, gof = 36)
# likelihood-ratio test
llike_rtest_pv <- pchisq(
q = 2*(ew_lrtn_m8$loglik-ew_lrtn_m7$loglik),
df = length(ew_lrtn_m8$coef)-length(ew_lrtn_m7$coef),
lower.tail = FALSE
)
# cat("Likelihood-Ratio test p-value:", llike_rtest_pv, "")
lrt_res <- lrtest(ew_lrtn_m7, ew_lrtn_m8)
sprintf(
"Likelihood-Ratio test p-value: %f (lrtest: %d; %f; %f)",
llike_rtest_pv,
lrt_res$Df[2],
lrt_res$Chisq[2],
lrt_res$`Pr(>Chisq)`[2]
)
In [745]:
ew_lrtn_m9 = arima(ew_lrtn_ts, order = c(0, 0, 1), seasonal = c(1, 0, 1), method = 'ML')
rsq = 1 - sum(ew_lrtn_m9$residuals^2)/sum(ew_lrtn^2)
summary(ew_lrtn_m9); rsq;
par(bg = 'white')
tsdiag1(ew_lrtn_m9, gof = 36)
# likelihood-ratio test
llike_rtest_pv <- pchisq(
q = 2*(ew_lrtn_m9$loglik-ew_lrtn_m7$loglik),
df = length(ew_lrtn_m9$coef)-length(ew_lrtn_m7$coef),
lower.tail = FALSE
)
# cat("Likelihood-Ratio test p-value:", llike_rtest_pv, "")
lrt_res <- lrtest(ew_lrtn_m7, ew_lrtn_m9)
sprintf(
"Likelihood-Ratio test p-value: %f (lrtest: %d; %f; %f)",
llike_rtest_pv,
lrt_res$Df[2],
lrt_res$Chisq[2],
lrt_res$`Pr(>Chisq)`[2]
)
2-13 (c)¶
In [714]:
npts = 12
h = 12
eotr = length(ew_lrtn)-h
freq = 12
order = c(1,0,0)
ar1_fc_res1 <- plot_arima_forecast_fig(
ew_lrtn_ts, eotr, h, npts, freq, gof=36,
order=order, method='ML'
)
ar1_fc_res2 <- plot_arima_forecast_fig(
ew_lrtn_ts, length(ew_lrtn), h, npts, freq, gof=36,
order=order, method='ML'
)
In [713]:
npts = 12
h = 12
eotr = length(ew_lrtn)-h
freq = 12
order = c(1,0,0)
seasonal = list(order = c(1,0,0), period = freq)
ar2_fc_res1 <- plot_arima_forecast_fig(
ew_lrtn_ts, eotr, h, npts, freq, gof=36,
order=order, seasonal=seasonal, method='ML'
)
ar2_fc_res2 <- plot_arima_forecast_fig(
ew_lrtn_ts, length(ew_lrtn), h, npts, freq, gof=36,
order=order, seasonal=seasonal, method='ML'
)
In [715]:
npts = 12
h = 12
eotr = length(ew_lrtn)-h
freq = 12
order = c(0,0,1)
ma1_fc_res1 <- plot_arima_forecast_fig(
ew_lrtn_ts, eotr, h, npts, freq, gof=36,
order=order, method='ML'
)
ma1_fc_res2 <- plot_arima_forecast_fig(
ew_lrtn_ts, length(ew_lrtn), h, npts, freq, gof=36,
order=order, method='ML'
)
In [716]:
npts = 12
h = 12
eotr = length(ew_lrtn)-h
freq = 12
order = c(0,0,1)
seasonal = list(order = c(0,0,1), period = freq)
ma2_fc_res1 <- plot_arima_forecast_fig(
ew_lrtn_ts, eotr, h, npts, freq, gof=36,
order=order, seasonal=seasonal, method='ML'
)
ma2_fc_res2 <- plot_arima_forecast_fig(
ew_lrtn_ts, length(ew_lrtn), h, npts, freq, gof=36,
order=order, seasonal=seasonal, method='ML'
)
In [726]:
npts = 12
h = 12
eotr = length(ew_lrtn)-h
freq = 12
order = c(0,0,1)
seasonal = list(order = c(1,0,1), period = freq)
arma1_fc_res1 <- plot_arima_forecast_fig(
ew_lrtn_ts, eotr, h, npts, freq, gof=36,
order=order, seasonal=seasonal, method='ML'
)
arma1_fc_res2 <- plot_arima_forecast_fig(
ew_lrtn_ts, length(ew_lrtn), h, npts, freq, gof=36,
order=order, seasonal=seasonal, method='ML'
)
2-13 (d)¶
- Likelihood ratio test.
In [741]:
print("1 vs 3")
lrtest(ew_lrtn_m1, ew_lrtn_m3);
print("2 vs 6")
lrtest(ew_lrtn_m2, ew_lrtn_m6);
print("3 vs 6")
lrtest(ew_lrtn_m3, ew_lrtn_m6);
print("3 vs 7")
lrtest(ew_lrtn_m3, ew_lrtn_m7);
print("6 vs 7")
lrtest(ew_lrtn_m6, ew_lrtn_m7);
print("7 vs 8")
lrtest(ew_lrtn_m7, ew_lrtn_m8);
print("1 vs 9")
lrtest(ew_lrtn_m1, ew_lrtn_m9);
print("3 vs 9")
lrtest(ew_lrtn_m3, ew_lrtn_m9);
print("7 vs 9")
lrtest(ew_lrtn_m7, ew_lrtn_m9);
print("8 vs 9")
lrtest(ew_lrtn_m8, ew_lrtn_m9);
In [ ]:
2-14 (S&P 500 Spot and Futures)¶
- This is a one-minute data. They are log-price.
- Several authors used the data to study index futures arbitrage.
- Interesting observations.
- The
lg(spot_rtn)has momentum inpacfandacf, butlg(fut_rtn)has no momentum inpacfandacf. This means the future is much more liquid than spot.
- The
In [451]:
da = read.table("../AFTS_sol/data/sp5may.dat", header = T)
dim(da); da[1:5, ]
In [452]:
lg_spt = da$lnspot
lg_fut = da$lnfuture
lg_spt_rtn = diff(lg_spt)
lg_fut_rtn = diff(lg_fut)
In [28]:
plot_time_fig(lg_spt, main = "Spot", xlab = "Minute")
plot_time_fig(lg_fut, main = "Futures", xlab = "Minute")
par(bg = 'white')
plot(lg_spt, type = 'l', col = 'blue', lty = 1)
lines(lg_fut, type = 'l', col = 'red', lty = 2)
legend("topleft", legend = c("Spot", "Futures"), col = c('blue', 'red'), lty = c(1, 2))
In [29]:
plot_time_fig(lg_spt_rtn, main = "Spot Return", xlab = "Minute")
plot_time_fig(lg_fut_rtn, main = "Futures Return", xlab = "Minute")
par(bg = 'white')
plot(lg_spt_rtn, type = 'l', col = 'blue', lty = 1)
lines(lg_fut_rtn, type = 'l', col = 'red', lty = 2)
legend("topleft", legend = c("Spot", "Futures"), col = c('blue', 'red'), lty = c(1, 2))
In [103]:
par(bg = 'white')
plot(lg_spt, lg_fut)
plot(lg_spt_rtn, lg_fut_rtn)
In [31]:
plot_pacf_acf(lg_spt)
In [32]:
plot_pacf_acf(lg_spt_rtn, main = 'Spot Rtn')
plot_pacf_acf(lg_fut_rtn, main = 'Fut Rtn')
Simple linear model without time series errors¶
In [37]:
library(lmtest)
In [87]:
snp_m0 = lm(lg_fut_rtn~1+lg_spt_rtn)
rsq = 1 - sum(snp_m0$residuals^2)/sum(lg_fut_rtn^2)
summary(snp_m0); rsq
In [88]:
snp_m1 = lm(lg_fut_rtn~-1+lg_spt_rtn)
rsq = 1 - sum(snp_m1$residuals^2)/sum(lg_fut_rtn^2)
summary(snp_m1); rsq
In [57]:
lrtest(snp_m0, snp_m1)
Arima model with time series errors¶
eacfshows(1, 0, 1)as thearimaorder.
In [84]:
library(forecast)
In [144]:
perform_and_print_eacf(snp_lm1$residuals, ar.max = 25, ma.max = 12)
In [40]:
cat("order =", ar(snp_lm1$residuals, method = 'mle')$order, "\n")
Use full data set¶
In [89]:
snp_m2 = arima(lg_fut_rtn, order = c(1, 0, 0), xreg = lg_spt_rtn, include.mean = F, method = "ML")
rsq = 1 - sum(snp_m2$residuals^2)/sum(lg_fut_rtn^2)
summary(snp_m2); rsq;
lrtest(snp_m1, snp_m2)
In [90]:
snp_m3 = arima(lg_fut_rtn, order = c(0, 0, 1), xreg = lg_spt_rtn, include.mean = F, method = "ML")
rsq = 1 - sum(snp_m3$residuals^2)/sum(lg_fut_rtn^2)
summary(snp_m3); rsq;
lrtest(snp_m1, snp_m3)
In [91]:
snp_m4 = arima(lg_fut_rtn, order = c(12, 0, 0), xreg = lg_spt_rtn, include.mean = F, method = "ML")
rsq = 1 - sum(snp_m4$residuals^2)/sum(lg_fut_rtn^2)
summary(snp_m4); rsq;
lrtest(snp_m2, snp_m4)
In [92]:
snp_m5 = arima(lg_fut_rtn, order = c(14, 0, 0), xreg = lg_spt_rtn, include.mean = F, method = "ML")
rsq = 1 - sum(snp_m5$residuals^2)/sum(lg_fut_rtn^2)
summary(snp_m5); rsq;
lrtest(snp_m4, snp_m5)
In [93]:
snp_m6 = arima(lg_fut_rtn, order = c(0, 0, 12), xreg = lg_spt_rtn, include.mean = F, method = "ML")
rsq = 1 - sum(snp_m6$residuals^2)/sum(lg_fut_rtn^2)
summary(snp_m6); rsq;
lrtest(snp_m4, snp_m6)
In [94]:
snp_m7 = arima(lg_fut_rtn, order = c(0, 0, 10), xreg = lg_spt_rtn, include.mean = F, method = "ML")
rsq = 1 - sum(snp_m7$residuals^2)/sum(lg_fut_rtn^2)
summary(snp_m7); rsq;
lrtest(snp_m4, snp_m7);
lrtest(snp_m6, snp_m7);
In [95]:
snp_m8 = arima(lg_fut_rtn, order = c(2, 0, 1), xreg = lg_spt_rtn, include.mean = F, method = "ML")
rsq = 1 - sum(snp_m8$residuals^2)/sum(lg_fut_rtn^2)
summary(snp_m8); rsq;
lrtest(snp_m6, snp_m8);
lrtest(snp_m7, snp_m8);
In [200]:
snp_m9 = auto.arima(
lg_fut_rtn,
d = 0,
D = 0,
max.p = 3,
max.q = 3,
max.P = 0,
max.Q = 0,
max.order = 6,
seasonal = TRUE,
method = "ML",
allowmean = TRUE,
xreg = as.matrix(lg_spt_rtn),
stepwise = FALSE,
parallel = TRUE,
num.cores = 4
)
rsq = 1 - sum(snp_m9$residuals^2)/sum(lg_fut_rtn^2)
summary(snp_m9); rsq;
lrtest(snp_m4, snp_m9);
lrtest(snp_m5, snp_m9);
lrtest(snp_m6, snp_m9);
lrtest(snp_m7, snp_m9);
lrtest(snp_m8, snp_m9);
In [142]:
snp_m9$loglik/length(lg_fut_rtn)
In [ ]:
Subsample data so that we focus more on larger magnitude of spot return¶
- Seems like sub-sampling near the timestamps that are large in magnitude gives higher R-squared, even though here the R-squared may not be a rigorious metric to use.
In [129]:
spt_qtl_90 = quantile(abs(lg_spt_rtn), probs = 0.9)
spt_qtl_90
In [139]:
byd_qtl_90_lg_spt_rtn_idx = which(abs(lg_spt_rtn)>spt_qtl_90)
byd_qtl_90_lg_spt_rtn_nei_idx = c()
for (i in 0:4) {
byd_qtl_90_lg_spt_rtn_nei_idx = c(byd_qtl_90_lg_spt_rtn_nei_idx, byd_qtl_90_lg_spt_rtn_idx + i)
}
print(length(byd_qtl_90_lg_spt_rtn_nei_idx))
byd_qtl_90_lg_spt_rtn_nei_idx = sort(unique(byd_qtl_90_lg_spt_rtn_nei_idx))
print(length(byd_qtl_90_lg_spt_rtn_nei_idx))
new_lg_fut_rtn = lg_fut_rtn[byd_qtl_90_lg_spt_rtn_nei_idx]
new_lg_spt_rtn = lg_spt_rtn[byd_qtl_90_lg_spt_rtn_nei_idx]
In [141]:
snp_m10 = auto.arima(
new_lg_fut_rtn,
d = 0,
D = 0,
max.p = 3,
max.q = 3,
max.P = 0,
max.Q = 0,
max.order = 6,
seasonal = TRUE,
method = "ML",
allowmean = TRUE,
xreg = as.matrix(new_lg_spt_rtn),
stepwise = FALSE,
parallel = TRUE,
num.cores = 4
)
rsq = 1 - sum(snp_m10$residuals^2)/sum(new_lg_fut_rtn^2)
summary(snp_m10); rsq;
In [143]:
snp_m10$loglik/length(new_lg_fut_rtn)
In [ ]:
Forecasting¶
- F1: Use the model fit with all data to forecast.
- F2: Use the model fit with subsampled data to forecast.
- The point-wise prediction are very close, b/c the coefs are close to each other.
- The F2 gives larger CI than F1, b/c the $\hat{\sigma}^2$ and all the variance of estimates of coefs are larger for F2 than F1. This is b/c in fitting F2, a data set with smaller sample size is used.
Use full data set¶
In [100]:
# auto.arima
npts = 20
eotr = length(lg_fut_rtn)-npts
h = npts
freq = 1
xreg = as.matrix(lg_spt_rtn)
snp_fc_res = plot_auto_arima_forecast_fig(
da_ts=lg_fut_rtn, eotr=eotr, h=h, npts=npts, frequency=freq,
xreg=xreg,
main="Forecasts from Regression with ARIMA(2,0,1) errors\nfor S&P 500",
xlab="Minute", ylab="S&P 500 Minute Data", ylim=c(-0.001, 0.001),
d = 0,
D = 0,
max.p = 3,
max.q = 3,
max.P = 0,
max.Q = 0,
max.order = 6,
seasonal = TRUE,
method = "ML",
allowmean = TRUE,
stepwise = FALSE,
parallel = TRUE,
num.cores = 12
)
In [101]:
snp_fc_tb = comb_forecast_res(snp_fc_res, lg_fut_rtn, eotr, freq)
snp_fc_tb
In [524]:
# forecast:Arima
npts = 20
eotr = length(lg_fut_rtn)-npts
h = npts
freq = 1
xreg = as.matrix(lg_spt_rtn)
order = c(2,0,1)
seas = NULL
normal_snp_fc_res = plot_arima_forecast_fig(
da_ts=lg_fut_rtn, eotr=eotr, h=h, npts=npts, frequency=freq,
xreg=xreg, gof=36,
main=NULL, #"Forecasts from ARIMA(4,2,5)\nfor GDP Deflator",
xlab="Minute", ylab="S&P 500 Minute Data", ylim=c(-0.001, 0.001),
order=order, seasonal=seas, include.mean=F, method='ML'
)
# summary(normal_snp_fc_res) # Too much output
In [517]:
normal_snp_fc_tb = comb_forecast_res(normal_snp_fc_res, lg_fut_rtn, eotr, freq)
normal_snp_fc_tb
In [532]:
# forecast:Arima
npts = 20
eotr = length(lg_fut_rtn)-npts
h = npts
freq = 1
xreg = as.matrix(lg_spt_rtn)
order = c(2,0,1)
seas = NULL
sp_sz = length(lg_fut_rtn)
xreg = as.matrix(lg_spt_rtn[1:(sp_sz-npts)])
normal_snp_fm = Arima(
lg_fut_rtn[1:(sp_sz-npts)],
order = order,
xreg = xreg, method = "ML", include.mean = F)
summary(normal_snp_fm)
fc_xreg = as.matrix(lg_spt_rtn[(sp_sz-npts+1):sp_sz])
forecast(normal_snp_fm, h = h, xreg = fc_xreg)
Subsample data so that we focus more on larger magnitude of spot return¶
In [271]:
attributes(snp_m10)
In [273]:
snp_m10_sub_spl = deep_copy(snp_fc_res$model)
snp_m10_sub_spl$coef = deep_copy(snp_m10$coef)
snp_m10_sub_spl$sigma2 = deep_copy(snp_m10$sigma2)
snp_m10_sub_spl$var.coef = deep_copy(snp_m10$var.coef)
In [287]:
display(snp_m10_sub_spl)
In [275]:
# Using the model fit with subsampled data
snp_sub_spl_ts_fc_res = forecast(
snp_m10_sub_spl,
h = h,
xreg = as.matrix(lg_spt_rtn[(length(lg_spt_rtn)-h+1):length(lg_spt_rtn)])
)
snp_sub_spl_ts_fc_res
In [284]:
# Using the model fit with all data
display(snp_fc_res)
In [280]:
# Forecast using the model fit with subsampled data
draw_arima_forecast_fig(
lg_fut_rtn, eotr=eotr, h=h, npts=npts, frequency=freq, ts_fc_res=snp_sub_spl_ts_fc_res,
main="Using model that fitted using more weights on timestamps\nwith 90-percentile magnitude spot return",
xlab=NULL, ylab=NULL, ylim = c(-0.001, 0.001)
)
In [286]:
# Forecast using the model fit with all data
draw_arima_forecast_fig(
lg_fut_rtn, eotr=eotr, h=h, npts=npts, frequency=freq, ts_fc_res=snp_fc_res,
main="Using model that fitted using more weights on timestamps\nwith all data (except the forecast part)",
xlab=NULL, ylab=NULL, ylim = c(-0.001, 0.001)
)
In [ ]:
2-15 (Quarterly GDP Implicit Prce Deflator)¶
- This is raw time series data.
- Seems that need to take diff twice to get a stationary time series.
In [317]:
require(xts, TSA, lmtest)
In [290]:
da = read.table("../AFTS_sol/data/q-gdpdef.txt", header = T)
dim(da); da[1:5,]
In [309]:
gdpdef = da$gdpdef
gdpdef_ts = ts(gdpdef, frequency = 4, start = c(1947, 1))
diff_gdpdef = diff(gdpdef)
diff_gdpdef_ts = ts(diff_gdpdef, frequency = 4, start = c(1947, 2))
lg_gdpdef = log(gdpdef)
lg_gdpdef_ts = ts(lg_gdpdef, frequency = 4, start = c(1947, 1))
diff_lg_gdpdef = diff(lg_gdpdef)
diff_lg_gdpdef_ts = ts(diff_lg_gdpdef, frequency = 4, start = c(1947, 2))
In [292]:
plot_time_fig(gdpdef_ts, main = "GDP Deflator", xlab = "Year")
plot_time_fig(lg_gdpdef_ts, main = "log(GDP Deflator)", xlab = "Year")
In [310]:
plot_time_fig(diff_gdpdef_ts, main = "diff(GDP Deflator)", xlab = "Year")
plot_time_fig(diff_lg_gdpdef_ts, main = "diff(log(GDP Deflator))", xlab = "Year")
In [311]:
plot_time_fig(diff(diff_gdpdef_ts), main = "diff(diff(GDP Deflator))", xlab = "Year")
plot_time_fig(diff(diff_lg_gdpdef_ts), main = "diff(diff(log(GDP Deflator)))", xlab = "Year")
In [305]:
lag.max = 24
plot_pacf_acf(gdpdef, lag.max = lag.max, xlim = c(0, lag.max))
plot_pacf_acf(lg_gdpdef, lag.max = lag.max, xlim = c(0, lag.max))
In [306]:
lag.max = 24
plot_pacf_acf(diff(gdpdef), lag.max = lag.max, xlim = c(0, lag.max))
plot_pacf_acf(diff(lg_gdpdef), lag.max = lag.max, xlim = c(0, lag.max))
In [ ]:
Fit models¶
In [541]:
perform_and_print_eacf(diff_lg_gdpdef, ar.max = 25, ma.max = 12)
In [314]:
perform_and_print_eacf(diff(diff_lg_gdpdef), ar.max = 25, ma.max = 12)
In [313]:
cat("order =", ar(diff(diff_lg_gdpdef), method = 'mle')$order, "\n")
In [319]:
gof = 24
In [320]:
gdpdef_m1 = arima(lg_gdpdef, order = c(4, 2, 3), method = "ML")
summary(gdpdef_m1)
tsdiag1(gdpdef_m1, gof.lag = gof)
In [322]:
gdpdef_m2 = arima(
lg_gdpdef,
fixed = c(NA, 0, NA, NA, NA, NA, NA),
order = c(4, 2, 3), method = "ML")
summary(gdpdef_m2)
tsdiag1(gdpdef_m2, gof.lag = gof);
lrtest(gdpdef_m1, gdpdef_m2)
In [323]:
gdpdef_m3 = arima(
lg_gdpdef,
order = c(6, 2, 5), method = "ML")
summary(gdpdef_m3)
tsdiag1(gdpdef_m3, gof.lag = gof);
lrtest(gdpdef_m1, gdpdef_m3)
In [325]:
gdpdef_fm <- auto.arima(
lg_gdpdef,
d = NA,
D = 0,
max.p = 6,
max.d = 3,
max.q = 6,
max.P = 0,
max.Q = 0,
max.order = 12,
seasonal = TRUE,
method = "ML",
allowmean = TRUE,
stepwise = FALSE,
parallel = TRUE,
num.cores = 12
)
summary(gdpdef_fm)
In [327]:
gdpdef_m4 = arima(
lg_gdpdef,
order = c(2, 2, 5), method = "ML")
summary(gdpdef_m4)
tsdiag1(gdpdef_m4, gof.lag = gof);
lrtest(gdpdef_m1, gdpdef_m4)
lrtest(gdpdef_m3, gdpdef_m4)
In [328]:
gdpdef_m5 = arima(
lg_gdpdef,
fixed = c(NA, NA, 0, NA, NA, 0, NA),
order = c(2, 2, 5), method = "ML")
summary(gdpdef_m5)
tsdiag1(gdpdef_m5, gof.lag = gof);
lrtest(gdpdef_m1, gdpdef_m5)
lrtest(gdpdef_m3, gdpdef_m5)
lrtest(gdpdef_m4, gdpdef_m5)
In [352]:
gdpdef_m6 = arima(
lg_gdpdef,
order = c(4, 2, 5), method = "ML")
summary(gdpdef_m6)
In [354]:
gdpdef_m7 = arima(
lg_gdpdef,
fixed = c(0, NA, 0, NA, NA, NA, NA, NA, NA),
order = c(4, 2, 5), method = "ML")
summary(gdpdef_m7)
tsdiag1(gdpdef_m5, gof.lag = gof);
lrtest(gdpdef_m5, gdpdef_m7)
In [ ]:
Forecast¶
In [573]:
npts = 8
eotr = length(lg_gdpdef_ts)-npts
h = npts
freq = 4
gdpdef_fc_res = plot_auto_arima_forecast_fig(
da_ts=lg_gdpdef_ts, eotr=eotr, h=h, npts=npts, frequency=freq,
main="Forecasts from ARIMA(4,2,5)\nfor GDP Deflator",
xlab="Year", ylab="lg(GDP Deflator)", ylim=c(4.7, 4.95),
d = NA,
D = 0,
max.p = 6,
max.d = 3,
max.q = 6,
max.P = 0,
max.Q = 0,
max.order = 12,
seasonal = TRUE,
method = "ML",
allowmean = TRUE,
stepwise = FALSE,
parallel = TRUE,
num.cores = 4
)
summary(gdpdef_fc_res)
In [634]:
gdpdef_fc_tb = comb_forecast_res(gdpdef_fc_res, lg_gdpdef_ts, eotr, freq)
gdpdef_fc_tb
In [550]:
gdpdef_fc_res1 = plot_auto_arima_forecast_fig(
da_ts=lg_gdpdef_ts, eotr=length(lg_gdpdef_ts), h=h, npts=npts, frequency=freq,
main="Forecasts from ARIMA(2,2,5)\nfor GDP Deflator",
xlab="Year", ylab="lg(GDP Deflator)", ylim=c(4.7, 4.95),
d = NA,
D = 0,
max.p = 6,
max.d = 3,
max.q = 6,
max.P = 0,
max.Q = 0,
max.order = 12,
seasonal = TRUE,
method = "ML",
allowmean = TRUE,
stepwise = FALSE,
parallel = TRUE,
num.cores = 4
)
summary(gdpdef_fc_res1)
In [ ]:
In [551]:
order = c(2,2,5)
normal_gdpdef_fc_res = plot_arima_forecast_fig(
da_ts=lg_gdpdef_ts, eotr=eotr, h=h, npts=npts, frequency=freq,
main=NULL, #"Forecasts from ARIMA(4,2,5)\nfor GDP Deflator",
xlab="Year", ylab="lg(GDP Deflator)", ylim=c(4.7, 4.95),
order=order, method='ML'
)
In [635]:
normal_gdpdef_fc_tb = comb_forecast_res(normal_gdpdef_fc_res, lg_gdpdef_ts, eotr, freq)
normal_gdpdef_fc_tb
In [585]:
order = c(2,2,5)
normal_gdpdef_fc_res1 = plot_arima_forecast_fig(
da_ts=lg_gdpdef_ts, eotr=length(lg_gdpdef_ts), h=h, npts=npts, frequency=freq,
main="Forecasts from ARIMA(2,2,5)\nfor GDP Deflator",
xlab="Year", ylab="lg(GDP Deflator)", ylim=c(4.7, 4.95),
order=order, method='ML'
)
In [616]:
comb_forecast_res(normal_gdpdef_fc_res1, lg_gdpdef_ts, length(lg_gdpdef_ts), freq)
In [589]:
gdpdef_fc_res$model$coef; gdpdef_fc_res$mean;
normal_gdpdef_fc_res$model$coef; normal_gdpdef_fc_res$mean
In [ ]:
Alternative solution¶
In [540]:
require(xts)
da = read.table("../AFTS_sol/data/q-gdpdef.txt", header = T)
y = log(da$gdpdef)
plot(xts(y, order.by = as.Date(paste(da$year, da$mom, da$day, sep = '-'))),
type = 'o', main = '', xlab = 'date', ylab = 'y')
par(mfrow = c(2, 1), bg = 'white')
pacf(diff(y))
acf(diff(y))
In [539]:
par(bg = 'white')
m1 = arima(diff(y), order = c(2, 0, 1))
m1
tsdiag(m1, gof = 20)
m2 = arima(diff(y), order = c(2, 0, 1), seasonal = list(order = c(0, 0, 1), period = 7))
m2
tsdiag(m2, gof = 20)
mp = predict(m2, 4)
Y1 = y[length(y)] + mp$pred[1]
Y2 = Y1 + mp$pred[2]
Y3 = Y2 + mp$pred[3]
Y4 = Y3 + mp$pred[4]
cat("prediction =", Y1, Y2, Y3, Y4, "\n")
In [569]:
npts = 8
h = npts
freq = 4
order = c(2,1,1)
seas = list(order = c(0,0,1), period = 7)
alt_normal_gdpdef_fc_res = plot_arima_forecast_fig(
da_ts=lg_gdpdef_ts, eotr=length(lg_gdpdef_ts), h=h, npts=npts, frequency=freq,
main=NULL, #"Forecasts from ARIMA(4,2,5)\nfor GDP Deflator",
xlab="Year", ylab="lg(GDP Deflator)", ylim=c(4.7, 4.95),
order=order, seasonal=seas, method='ML', include.drift=T
)
In [571]:
alt_normal_gdpdef_fc_res$model;
m2;
alt_normal_gdpdef_fc_res$mean
In [572]:
cat("prediction =", Y1, Y2, Y3, Y4, "\n")
In [ ]:
End¶
In [ ]:
In [ ]:
In [ ]: